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We present a statistical equilibrium model of self-organization in a class of focusing, nonintegrable 
nonlinear Schrodinger (NLS) equations. The theory predicts that the asymptotic-time behavior 
of the NLS system is characterized by the formation and persistence of a large-scale coherent 
solitary wave, which minimizes the Hamiltonian given the conserved particle number (L'^-norm 
squared), coupled with small-scale random fluctuations, or radiation. The fluctuations account 
for the difference between the conserved value of the Hamiltonian and the Hamiltonian of the 
coherent state. The predictions of the statistical theory are tested against the results of direct 
numerical simulations of NLS, and excellent qualitative and quantitative agreement is demonstrated. 
In addition, a careful inspection of the numerical simulations reveals interesting features of the 
transitory dynamics leading up to the to the long-time statistical equilibrium state starting from a 
given initial condition. As time increases, the system investigates smaller and smaller scales, and it 
appears that at a given intermediate time after the coalescense of the soliton structures has ended,, 
the system is nearly in statistical equilibrium over the modes that it has investigated up to that 
time. 
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I. INTRODUCTION: NLS AND SOLITON TURBULENCE 



A fascinating feature of many turbulent fluid and plasma systems is the emergence and persistence of large-scale 
organized states, or coherent structures, in the midst of small-scale turbulent fluctuations. A familiar example is 
the formation of macroscopic quasi-steady vortices in a turbulent large Reynolds number two dimensional fluid [0-^ . 
Such phenomena also occur for many classical Hamiltonian systems, even though the dynamics of these systems 
is formally reversible In the present work, we shall focus our attention on another class of nonlinear partial 
differential equations whose solutions exhibit the tendency to form persistent coherent structures immersed in a sea of 
microscopic turbulent fluctuations. This is the class of nonlinear wave systems described by the well-known nonlinear 
Schrodinger (NLS) equation: 

+ AV' + /(l^nv- =0, (1) 

where ip{r,t) is a complex field and A is the Laplacian operator. The NLS equation describes the slowly- varying 
envelope of a wave train in a dispersive conservative system. It models, among other things, gravity waves on deep 
water ||^, Langmuir waves in plasmas j^], pulse propagation along optical fibers |Q, and superfluid dynamics 
When /(IV'P) = 8'iid eqn. (|l|) is posed on the whole real line or on a bounded interval with periodic boundary 

conditions, the equation is completely integrable Otherwise, it is nonintegrable. 

The NLS equation (|l|) may be cast in the Hamiltonian form idtijj = SH/6ip*, where ip* is the complex conjugate of 
the field ip, and H is the Hamiltonian: 

H{^)= J (iVVp-F(IVn) dr. (2) 

Here, the potential F is defined via the relation F{a) — f{y) dy. The dynamics (^) conserves, in addition to the 
Hamiltonian, the particle number 

N{i;)= [ \^\^dr. (3) 



We shall assume throughout that eqn. (|^) is posed in a bounded one dimensional interval with either periodic 
or homogeneous Dirichlet boundary conditions. We restrict our attention to attractive, or focusing, nonlinearities 
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/(/(a) > 0, /'(a) > 0) such that the dynamics described by (||) is nonintegrable, free of wave collapse, and admits 
stable solitary-wave solutions. The dynamics under these conditions has been referred to as soliton turbulence [0. 
Such is the case for the important power law nonlinearities, /(IV'P) = iV'l^j with < s < 4 (in the periodic case, s ^ 2 
for nonintegrability) and also for the physically relevant saturated nonlinearities = + IV'P) ^^^d 

/(I^P) = 1 ~ exp(— l^jp), which arise as corrections to the cubic nonlincarity for large wave amplitudes JlBf . 

Equation in one spatial dimension has solitary wave solutions of the form ip{x,t) — (j){x) exp{iX^t), where <j) 
satisfies the nonlinear eigenvalue equation: 



0. 



(4) 



It has been argued ||T^,0 that the solitary wave solutions play a prominent role in the long-time dynamics of (|l|) , in 
that they act as statistical attractors to which the system relaxes. The numerical simulations in [lO[ |, as well as the 
simulations we shall present within this article, support this conclusion. Indeed, it is seen that for rather generic initial 
conditions the field ip evolves, after a sufficiently long time, into a state consisting of a spatially localized coherent 
structure, which compares quite favorably to a solution of (^, immersed in a sea of turbulent small-scale turbulent 
fluctuations. At intermediate times the solution typically consists of a collection of these soliton-like structures, but 
as time evolves, the solitons undergo a succession of collisions in which the smaller soliton decreases in amplitude, 
while the larger one increases in amplitude. When solitons collide or interact, they shed radiation, or small-scale 
fluctuations. The interaction of the solitons continues until eventually a single soliton of large amplitude survives 
amidst the turbulent background radiation. Figure (|^) below illustrates the evolution of the solution of (|l|) for the 
particular nonlinearity /(IV'P) — IV'I and with periodic boundary conditions on the spatial interval [0, 256]. 
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FIG. 1. Profile of the modulus It/j^ at four different times for the system (0) with nonlinearity /(IV'P) ~ IV'I periodic 
boundary conditions on the interval [0,256]. The initial condition is ^{x,t = 0) = A, with A — 0.5, plus a small random 
perturbation. The numerical scheme used to approximate the solution is the split-step Fourier method. The grid size is 
dx = 0.125, and the number of modes is n = 2048. a) t = 50 unit time: Due to the modulational instability, an array of 
soliton-like structures separated by the typical distance U = 2n / A/2 = 4-k is created; h) t = 1050 unit time: The solitons 
interact and coalesce, giving rise to a smaller number of solitons of larger amplitude; c) t — 15050: The coarsening process has 
ended. One large soliton remains in a background of small-amplitude radiation. Notice that for t — 55050 unit time (Figure 
d)), the amplitude of the fiuctuations has diminished while the amplitude of the soliton has increased. 

In modeling the long-time behavior of a Hamiltonian system such as NLS, it seems natural to appeal to the methods 
of equilibrium statistical mechanics. That such an approach may be relevant for understanding the asymptotic-time 
state for NLS has already been suggested in ]|l^ , although the thermodynamic arguments presented by these authors 
are rather formal and somewhat incomplete. Motivated in part by the ideas outlined in po| , Jordan et al. psf ] have 
recently constructed a mean-field statistical theory to characterize the large-scale structure and the statistics of the 
small-scale fluctuations inherent in the asymptotic-time state of the focusing nonintegrable NLS system (fl). The 
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main prediction of this theory is that the coherent state that emerges in the long-time hmit is the ground state 
solution of equation (^. That is, it is the solitary wave that minimizes the Hamiltonian H given the constraint 
N — N''\ where N'^ is the initial and conserved value of the particle number integral. This prediction is in accord 
with previous theories |lcHl^ , but the approach taken in |l5[] is new, and provides a definite interpretation to the 
notion set forth in the earlier works that it is "thermodynamically advantageous" for the NLS system to approach 
a coherent solitary wave structure that minimizes the Hamiltonian subject to fixed particle number. The statistical 
theory also gives predictions for the particle number spectral density and the kinetic energy spectral density, at least 
for a finite-dimensional spectral truncation of the NLS dynamics (|l|). In particular, it predicts an equipartition of 
kinetic energy among the small-scale fluctuations. 

In the present work, we shall begin with a brief review of this statistical theory. The predictions of the statistical 
theory will then be compared in detail with the results of direct numerical simulations of the NLS system. In addition, 
we will also closely examine the evolution of the particle number spectrum in our numerical simulations, as well as 
the dynamics (of finite spectral approximations) of the integrals Smi'ip) = J \D"^ip\'^ dx (here D™ denotes the m-th 
derivative with respect to the spatial variable). The statistical model, being strictly an equilibrium theory, does not 
give predictions concerning the finite time dynamics of these quantities. However, we shall see that it does give accurate 
estimates for the long-time saturation values of these quantities for a finite dimensional spectral approximation of 
the NLS dynamics. In addition, we will demonstrate that the integrals Sm exhibit interesting power law growth in 
time, as suggested by the weak turbulence theory developed by Pomeau jl^. 



II. MEAN-FIELD STATISTICAL MODEL 



In order to develop a meaningful statistical theory, we begin by introducing a finite-dimensional approximation 
of the NLS equation ([|). To fix ideas and notation, we will consider the NLS system with homogeneous Dirichlet 
boundary conditions on an interval fl of length L. Our methods can easily be modified to accommodate other 
boundary conditions, and we will consider below the predictions of the theory for periodic boundary conditions, as 
well. In addition, our techniques can easily be extended to higher dimensions, but we wish to concentrate on the 
one-dimensional case for ease of presentation. 

Let ej{x) = \/2/L siii{kjx) with kj — ttj'/L, and for any function g{x) on f2 denote by gj = J^g{x)ej{x) dx 
its jth Fourier coefficient with respect to the orthonormal basis ej,j = 1,2,- • •. Define the functions w*^"^(a;,i) = 
Xlj=i ^j(^)ei(a^) and v^^^x, t) — Y^^=i where the real coefficients Uj,Vj,j — 1, • • • , n, satisfy the coupled 

system of ordinary differential equations 

ii, - k^Vj + (/((u("))2 + (t;("))2)t,("))^. = 

(5) 

ilj + k^Uj - {jiiu^'^'^f + (f ("))2)u("))^. = . 

Then the complex function -0^") = u^") + iti(") satisfies the equation 

zV;^")+4",)+P"(/(|V^(")p)V'(")) = 0, 

where P" is the projection onto the span of the eigenfunctions ei, • • • ,e„. This equation is a natural spectral ap- 
proximation of the NLS equation ([^), and it may be shown that its solutions converge as n ^ cxd to solutions of (|l|) 

For given n, the system of equations (g) defines a dynamics on the 2n-dimensional phase space R^". This finite- 
dimensional dynamical system is a Hamiltonian system, with conjugate variables Uj and Wj, and with Hamiltonian 

= x„ + e„ , (6) 

where 




is the kinetic energy, and 

e„ = -i^F((ii("))' + («("))')dx, (8) 
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is the potential energy. The Hamihonian _ff„ is, of course, an invariant of the dynamics. The truncated version of the 
particle number 




is also conserved by the dynamics (jSV. The factor 1/2 is included in the definition of the particle number for 
convenience. The Hamiltonian system satisfies the Liouville property, which is to say that the measure IljLi dujdvj 
is invariant under the dynamics p7[ | . This property together with the assumption of ergodicity of the dynamics provide 
the usual starting point for a statistical treatment of a Hamiltonian system |l8| . 

With the finite dimensional Hamiltonian system in hand, we now consider a macroscopic description in terms 
of a probability density • • • , u„, vi • • • , Vn) on the 2n-dimensional phase-space R^". We seek a probability 

density that describes the statistical equilibrium state for the truncated dynamics. In accord with standard statistical 
mechanics and information theoretic principles, we define this state to be the density p^"^ on 2n-dimensional phase 
space which maximizes the Gibbs-Boltzmann entropy functional 

» n 

S{p) =- plogp Jl dujdvj , (10) 
subject to constraints dictated by the conservation of the Hamiltonian and the particle number under the dynamics 

(i ES- 

The usual canonical ensemble 



p (X exp {-PHn - pNn) , 

results from maximizing the entropy subject to the mean constraints (i?„) — and (iV„) — , where and are 
the given values of the Hamiltonian and the particle number, respectively, and [3 and /i are the Lagrange multipliers to 
enforce these constraints. However, it has been shown in ||lq , |20| that, for the focusing nonlinearities we consider here, 
the canonical ensemble is ill-defined in the sense that it is not normalizable (i.e., /p^2„ exp[— — pNn] 0^=1 dujdvj 
diverges ). Thus, wc are obliged to consider an alternative statistical equilibrium description of the NLS system based 
on constraints other than those that give rise to the canonical ensemble. The key to constructing an appropriate 
statistical model is based on the observation from numerical simulations that, for a large number of modes n, in 
the long-time limit, the field (u'^"\t;(")) decomposes into two essentially distinct components: a large-scale coherent 
structure, and small-scale radiation, or fluctuations. As time progresses, the amplitude of the fluctuations decreases, 
until eventually the contribution of the fluctuations to the particle number and the potential energy component of 
the Hamiltonian becomes negligible compared to the contribution from the coherent state, so that Nn and 8„ are 
determined almost entirely by the coherent structure. We have checked that this effect becomes even more pronounced 
when the resolution of the numerical simulations is improved (i.e., when the number of modes is increased with the 
length L of the spatial interval fixed). On the other hand, as the fluctuations exhibit rapid spatial variations, the 
amplitude of their gradient does not, in general, become negligible in the asymptotic time limit. Consequently, the 
fluctuations can make a signiflcant contribution to the kinetic energy component Kn of the Hamiltonian. This is 
illustrated in Fig. (||). 
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FIG. 2. Numerical simulation for the saturated non-linearity /(I'f/'l^) = + and for periodic boundary conditions. 

The total number of modes is n = 1024 and the spatial grid size is dx = 0.25, so that the length of periodic interval is L = 256. 
Displayed are the modulus of the field \ip\'^ (first and second rows), and the modulus of the gradient of the field \ipxf (third 
row) at unit times t = 30,000 (left) and t — 220,000 (right). The second row shows the same results as the first row, except 
that the we have restricted the range on the vertical axis in order to focus in on the the fluctuations of the field. Notice that 
the dynamics for this saturated nonlinearity is qualitatively similar to that for the power law nonlinearity /(IV'P) ~ \'^\ shown 
in Fig. (0): the long-time state consists of large-scale coherent solitary wave-like structure interacting with a sea of small-scale 
fluctuations (top row). The typical amplitude of the fluctuations of the field has decreased from t = 30,000 to t — 220,000 
(second row), while the amplitude of the coherent structure has increased somewhat. The maximum of the modulus of the field 
is on the order of 50 times larger than the typical modulus of the fluctuations a.t t — 220, 000. On the other hand, the typical 
amplitude of the fluctuations of the gradient of the held has actually increased somewhat from t — 30, 000 to t = 220, 000, 
and the typical amplitude of the fluctuations of the gradient is only several times smaller than the maximum amplitude of 
the gradient of the field (bottom row). Clearly, the fluctuations make a signiflcant contribution to the kinetic energy in the 
long-time limit. 



Denoting by (uj) and (vj) the means of the variables uj and vj with respect to the yet to be determined ensemble p^"-* , 
we now identify the coherent state with the mean-field pair ((^("'(a;)), (w*^"^(a;))) = {'^^^i{uj)ej{x),'^^^^{vj)ej{x)). 



The fluctuations, or 



small-scale radiation inherent in the 



between the state vector w^" 



long-time state then correspond to the difference 
) and the mean-field vector. The 



statistics of the fluctuations are encoded in the probability density p*^"' . Based on the considerations of the preceding 
paragraph, and the results of the numerical simulations displayed in Fig. (^, it seems reasonable to conjecture that 
the amplitude of the fluctuations of the field V'*-"'' in the long-time state of the NLS system (||) should vanish entirely 
(in some appropriate sense) in the continuum limit n — )■ oo. Thus we are led to the following vanishing of fluctuations 
hypothesis: 



(((5u("))2) + (((5w("))2) 



dx = 



[{iSu,r) + {{Sv,r)]^0, as 



(11) 



Here, Suj — Uj — (uj) represents the fluctuations of the Fourier coefficient Uj about its mean value (uj), and similarly 
for Svj. We emphasize that ( |ll|) is a hypothesis used to construct our statistical theory, and not a conclusion drawn 
from the theory itself. 

An immediate consequence of the vanishing of fluctuations hypothesis is that for n sufficiently large, the expectation 
(Nn) of the particle number is determined almost entirely by the mean ((u*^")), (tj^"))). Furthermore, the hypoth- 



ec ((« 



esis ( pT| ) implies that for n large, the expectation (e„(u 



(n) y(n) 



')) of the potential energy is well approximated by 



)), which is the potential energy of the mean. This may be seen by expanding the potential F about 



the mean ((u'^"^), (w^"^)) in equation (^, taking expectations, and noting that because of the vanishing of fluctuations 
hypothesis (plj), there holds |(e„(M("\ w^"))) — 6„((m(")), = o(l) as n ^ oo. Notice, however, that the vanish- 

ing of fluctuations hypothesis does not imply that the contribution of the fluctuations to the expectation of the kinetic 
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energy becomes negligible in the limit n — > oo. Indeed, this contribution is (1/2) + ii^'^j)'^)]' which 

need not tend to as n — > oo, even if (|ll|) holds. Thus, from these arguments, we conclude that for n sufficiently large, 
(Hn) « + (^i)) ~ h + (w^"^)^)dx. These considerations lead us to impose the following 

mean-field constraints on the admissible probability densities p on the 27i-dimensional phase space: 

^«(p)^iE;=i((".>' + (^.>') = ^^° 

(12) 

Hnip) ^ kT.U + («|)) -U^F{{u(-)f + {v^-))^)dx = . 

Here, and H'^ are the conserved values of the particle number and the Hamiltonian, as determined from initial 
conditions. The statistical equilibrium states are then defined to be probability densities p'-"-' on the phase-space R^" 
that maximize the entropy ( p^ ) subject to the constraints (|l^). We shall refer to the constrained maximum entropy 
principle that determines the statistical equilibria as (MEP)~^ 

Further justification and motivation for the vanishing of fluctuations hypothesis (|ll|), which leads to the mean-field 
constraints in the maximum entropy principle (MEP), are provided in [|l5|. In particular, it is proved in that 
the solutions p^"' of (MEP) concentrate on the phase-space manifold on which = and 7V„ A^^n the 
continuum limit n — > oo, in the sense that {Nn) — > A^", (iJ„) and var Nn 0, var iJ„ ^ in this limit. Here, 

var W denotes the variance of the random variable W. This concentration property establishes a form of asymptotic 
equivalence between the mean-field ensembles p'-"-' and the microcanonical ensemble, which is the invariant measure 
concentrated on the phase-space manifold on which iJ„ = iJ" and Nn = A^". It therefore provides a strong theoretical 
justification for the mean-field statistical model. 



III. CALCULATION AND ANALYSIS OF EQUILIBRIUM STATES 

The solutions p^"-' of (MEP) are calculated by an application of the Lagrange multiplier rule 

5'(p(")) = M<(p^"^)+/3<(p^"^), 

where /3 and /i are the Lagrange multipliers to enforce that the probability density p'^"^ satisfy the constraints (|l3). 
A straightforward but tedious calculation yields the following expression for the maximum entropy distribution p'^^ 



where, for j = 1, . . . , n. 



p^"\ui,...,Un,Vi,...,Vn) = J| (Uj , ) , (13) 



Pjiuj,vj) = ^ exp <^ ^ {{u, ~ {u,)f + {v, - {v,)f) \ , (14) 



with: 



(15) 

- ^ (/((^("))^ + («("))^)(t.(")))^. - ^(«,) . 

Thus, for each j, Uj and Vj are independent Gaussian variables, with means given by the nonlinear equations ( [T5| ) 
and with identical variances 

var Uj = var Vj — — . (16) 
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Note that var Uj = {{SujY) by definition, and likewise for Vj. Obviously, the multiplier (3 must be positive. Notice also 
that, since the probability density p'"' factors according to (p^), the Fourier modes u^, Vj^j = 1, • • • , n, are mutually 
uncorrelated. In addition, we see from ( p^ ) that the complex mean- field (V'^"^) = is solution of (setting 

A = M//3) 
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(V'^"^).. + P" - A(V>)) - 0, (17) 

which is clearly the spectral truncation of the eigenvalue equation for the continuous NLS system (|l|). It follows, 
therefore, that the mean-field predicted by our theory corresponds to a solitary wave solution of the NLS equation. 
Alternatively, the mean ((m^"-'), (w'^"^)) is a solution of the variational equation 5Hn + XSNn = 0, where A is a Lagrange 
multiplier to enforce the particle number constraint Nn = iV^. 

Now, as the maximum entropy distribution p*-"-* is required to satisfy the mean-field Hamiltonian constraint (|l^), 
it follows from (|l|)-(|l3) that 

i?o = ^ + H„ ((it (")),(!;("))). (18) 

The term n//3 represents the contribution to the kinetic energy from the Gaussian fluctuations, and i/„((u("'), (?;("))) 
is the Hamiltonian of the mean. Notice that the contribution of the fluctuations to the kinetic energy is divided evenly 
among the n Fourier modes. From dlq), we obtain the following expression for f3 in terms of the number of modes n 
and the Hamiltonian i/„((M'^"^), (w'^'T) of the mean: 

ifO-i/„((w(") ),(«("))) ■ ^^^^ 

Using equations (^3|)-([l^), we may easily calculate the entropy of any solution p*^"-* of (MEP). This yields, after 
some algebraic manipulations, that 



S(.<">) = C(„) + nlog ("''^°-^"f°'>-<°'"'»') 



(20) 



where C{n) = n — "^^^^ log(j^7r/2) depends only on the number of Fourier modes n. Clearly, the entropy S'(p^"^) will 
be maximum if and only if the mean field pair ((u*^"-*), (w^"))) corresponding to /d^"-* realizes the minimum possible 
value of Hn over all fields (u("\ ?;(")) that satisfy the constraint iV„('u("), -y^")) = N°. 

Equation ( po|) reveals that in statistical equilibrium the entropy is, up to additive and multiplicative constants, 
the logarithm of the the kinetic energy contained in the turbulent fluctuations about the mean state. This result, 
therefore, provides a precise interpretation to the notions set forth by Zakharov et al. ||l^ and Pomeau |14| that the 



entropy of the NLS system is directly related to the amount of kinetic energy contained in the small-scale fluctuations, 
and that it is "thermodynamically advantageous" for the solution of NLS to approach a ground state which minimizes 
the Hamiltonian for the given number of particles. 

We now know that iJ„((w(")), (w^"))) = H*, where H* is the minimum vale of i7„ allowed by the particle number 
constraint = As a consequence, the Lagrange multiplier /3 is uniquely determined by (p^): 

That the "inverse temperature" (3 scales linearly with the number of Fourier modes n is required in order to obtain a 
meaningful continuum limit n — > oo in which the expectations of the Hamiltonian and particle number remain finite. 
The scaling of the inverse temperature with the number of modes is a common feature of the equilibrium statistical 
mechanics of finite dimensional approximations of other plasma and fluid systems with infinitely many degrees of 
freedom, as well pl| . The parameter A (which depends on n) is also determined by the requirement that the mean 
((«(")),(«("))) realize the minimum value of the Hamiltonian i7„ given the particle number constraint Nn = N^. 

Using eqns. ( |l6|) and (^l|), we may now obtain an exact expression for the contribution of the fluctuations to the 
expectation of the particle number. This is 

1 rrO rr* "'1 1 

-J2[{(Su,r) + {iSv,r)]^^-^J2p-Oi-), asn-oo. (22) 



Recall that in the derivation of the mean-field constraints (12|) , we assumed the vanishing of fluctuations condition 



). The calculation (P^ ) shows, therefore, that the maximum entropy distributions p^"-* indeed satisfy the hypothesis 
), and hence, that the mean-field statistical theory is consistent with the assumption that was made to derive it. 
But as the analysis of this section has shown, the maximum entropy distributions p^"-* provide much more information 
than is contained in the hypothesis (pi]). Most importantly, we know that the mean-field corresponding to p^"^ is 
an absolute minimizer of the Hamiltonian iJ„ subject to the particle number constraint Nn — N^. In addition, the 
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theory yields predictions for the particle number and kinetic energy spectral densities, at least for the 2n-dimensional 
spectrally truncated NLS system (H) with n large. Indeed, we have the following prediction for the particle number 
spectral density 

m') = m? + ^^^^^^ (23) 

where we have used the identity — Uj + ivj, and eqns. and The first ter m on the right hand side of ( |2^ ) 
is the contribution to the particle number spectrum from the mean, and the second term is the contribution from 
the fluctuations. Since the mean field is a smooth solution of the ground-state equation, its spectrum decays rapidly, 
so that for j >> 1, we have the approximation ~ — H*)/{nk'j). The kinetic energy spectral density is 

obtained simply by multiplying eqn. (j2^) by fc^. As emphasized above, we have the prediction that the kinetic energy 
arising from the fluctuations is equipartitioned among the n spectral modes, with each mode contributing the amount 
iH°-H*)/n. 

While we have chosen to present the statistical theory specifically for homogeneous Dirichlet boundary conditions, 
it is straightforward to develop the theory for NLS on a periodic interval of length L, as well. In this case, it is most 
convenient to write the spectrally truncated complex field ■0'^"^ as 

n/2 
j=-n/2 

for n an even positive integer, where kj ~ 2Trj / L. The predictions of the statistical theory remain the same as in 
the case of Dirichlet boundary conditions. In particular, the mean field (V-^")) is a minimizer of the Hamiltonian 
Hn given the particle number constraint iV„ = A^", and the particle number spectrum satisfies ( p3| ) for j ^ The 
Fourier coefficient ipQ may be consistently chosen to be deterministic (i.e., var -00 = and ("00) = 0o)i to eliminate 
the ambiguity arising from the mode.. 



IV. NUMERICAL RESULTS 



The general predictions of the statistical theory developed above do not depend crucially on the particular nonlin- 
earity / in the NLS equation (|l|). Indeed, for any / satisfying the conditions stated in the introduction, the coherent 
structure predicted by the theory in the continuum limit n —^ oo corresponds to the solitary wave that minimizes 
the energy for the given number of particles A^*^. Also, for any such nonlinearity /, the particle number spectrum in 
the long-time limit for the spectrally truncated NLS system (^, according to the statistical theory, should obey the 
relation (|2^). Of course, the minimum value H* of the Hamiltonian _ff„ which enters this formula does depend on /. 

Here, we will present numerical results primarily for periodic boundary conditions and for the focusing power law 
nonlinearity /(|0p) = |0|. That is, we shall solve numerically the particular NLS equation 

z9t0 + 9,,0+ 1010-0, (24) 

on a periodic interval of length L. We have, however, carried out similar numerical experiments for different focusing 
nonlinearities and for Dirichlet boundary conditions, and we observed that the general qualitative features of the long- 
time dynamics are unaltered by such changes. The nonlinearity /(|0p) = |0| actually represents a nice compromise 
between the focusing effect and nonlinear interactions. For weaker nonlinearities (such as the saturated ones), the 
interaction between modes is weak, and the time required to approach an asymptotic equilibrium state is quite long. 
On the other hand, for stronger nonlinearities, the solitary wave structures that emerge exhibit narrow peaks of large 
amplitude, and therefore, greater spatial resolution is required in the numerical simulations. 

The numerical scheme that we use for solving ( p^ is the well-known split-step Fourier method for a given number 
n of Fourier modes. Throughout the duration of the simulations, the relative error in the particle number is kept 
at less than 10~^ percent, and the relative error in the Hamiltonian is no greater than 0.1 percent. Notice that the 
numerical simulations, performed naturally for a finite number of modes, provide an ideal context for comparisons 
with the mean-field statistical theory outlined above. 

On the whole real line, the nonlinear Schrodinger equation ( p4[ ) has solitary wave solutions of the form 0(a::,<) = 
0(a;)e*-^ *, with 



2cosh^( 



(25) 
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The particle number N and the Hamiltonian H of these sohton-hke solutions are determined by the parameter A 
through the relationships N = 6A'^ and H = — -^A^. These solutions are centered at x = Xq, as shown in Figure (||), 
and because of the focusing property of equation (|2j), as N increases, the amplitude of the solitary wave increases, 
while its width decreases. For a given value of the particle number N, the solitary wave (|25| ) is the global minimizer 
of the Hamiltonian H (when the integrals in the definitions and ^ of the Hamiltonian and the particle number 
extend over the real line). Of course, the solitary wave solutions for the equation (|2j) on a finite interval, as well as 
those for the spectrally-truncated version (||), differ from the solution ( psf ) over the infinite interval. However, because 
the solitary waves (|2g) exhibit an exponential decay, for a large enough interval, and for a large enough number of 
modes n, such differences can be neglected for all practical purposes. 
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FIG. 3. Solitary wave solutions of equation ( ^ ) for particle numbers A*' = 1, 5 and 10. The modulus \tp\^ is shown as a 
function of the position. 

For constant A, the condensate il^{x,t) = Ae'^^* is an equilibrium solution of (p4[). However, since the nonlinearity is 
focusing, this spatially homogeneous solution is modulationally unstable. Indeed, if we expand tp around this solution 
in a scries of the form 

we obtain the dispersion relation: 

= Ae - 

Thus, the condensate is stable for > A, and unstable for k^ < A. The most unstable wave-number is ki = A/2. 

We choose to present in this paper the following set of numerical simulations: starting with the spatially ho- 
mogeneous solution i}}[x,t = 0) = A (with A of order 1), we add initially a small spatially uncorrelated random 
perturbation, so that the modulational instability develops. Although wc have checked that the long-time behavior 
of the solution is not dependent on the initial conditions, except through the initial and conserved values N'^ and 
of the particle number and the Hamiltonian, this class of initial conditions is particularly convenient for our purposes. 
For example, by considering different realizations of the initial random perturbation, we may perform an ensemble 
average over different initial conditions for a given A (and therefore for fixed and H^). Such initial conditions 
provide interesting analogies to standard fluid turbulence problems, as we will emphasize in the conclusion. 

The spatially uniform initial conditions we consider here may be thought of as being far away from the expected 
statistical attractor described by the maximum entropy probability density p^^\ Indeed, the spectrum of the con- 
densate differs considerably from the predicted statistical equilibrium spectrum (p3|). The numerical simulations that 
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we perform here provide strong evidence that the solutions of the spectrally truncated NLS system converge in the 
long-time limit to a state that may be considered as statistically steady. We shall compare the statistical properties of 
this long-time state with the predictions of the mean-field statistical theory that was developed and analyzed above. 
In addition, we shall also investigate the following questions concerning the nature of the evolution leading from the 
initial state to the long-time statistical equilibrium state: 1- How long does it take for the system to reach the vicinity 
of its statistical attractor, so that subsequently its statistical features may be considered as stationary? 2- How well 
can we characterize the "path" that a solution follows en route to the statistically steady state? That is, what are 
the generic features of the transitory dynamics? 

Figure (|l|) demonstrates that the transitory dynamics can be roughly decomposed into three stages: in the first 
stage, illustrated in Figure |l| a) , the modulational instability creates an array of soliton-like structures separated by a 
typical distance li — 2t: /ki associated with most unstable wave number ki. The second stage is characterized by the 
interaction and coalescence of these solitons. In this coarsening process, the number of solitons decreases, while the 
amplitudes of the surviving solitons increase, until eventually a single soliton of large amplitude persists amongst a sea 
of small-amplitude background radiation (Figures |]b) and c)). This intermediate stage has previously been observed 
for other nonlinear Schrodinger equations in one and two spatial dimensions |p^ , p2| , and it was shown in ||2^ that 
this coarsening process follows a self-similar dynamics. The dynamical exponents of these processes are not very well 
understood at this point, however. During the final stage of the dynamics, the surviving large-scale soliton interacts 
with the small-scale fluctuations. As time increases, the amplitude of the soliton increases, while the amplitude of 
the fluctuations decreases (note the changes from Figure |l| c) to Figure d)). In this stage of the dynamics, the mass 
(or number of particles) is gradually transferred from the small-scale fluctuations to the large-scale coherent soliton. 
For a flnite number of modes n, the dynamics eventually reaches a "stationary" state whose properties are very well 
described by the mean-fleld statistical equilibrium theory developed above, as we shall demonstrate. This implies that 
long-time state may, in fact, be thought of as a "statistical attractor", in the sense that, according to the statistical 
theory, it corresponds to a maximizer of the entropy functional ( p^ ) subject to the dynamical constraints (|l^). Note 
that because the dynamics is reversible, intermediate states such as those in Figure (|l] b) theoretically could still be 
attained even after the statistical equilibrium state has been reached. In fact, a numerical simulation starting from 
the state in Fig. (^ d)) but with the time step taken negative shows the reverse dynamics up to round-off errors, 
where one can observe the decomposition of the solution into an array of soliton-like structures as in figure (0 a)) 
for intermediate times, while in the limit t — oo an equilibrium state such as the one of figure (|l| d)) is once again 
attained. 

The tendency of the solution of the NLS system ( p4[ ) to approach the statistical equilibrium state is also captured in 
the evolution of the kinetic and potential energies (see Figure While the sum of these two quantities, which is the 
Hamiltonian, remains constant in time, we observe that the kinetic energy increases monotonically, and, consequently, 
the potential energy decreases monotonically as time goes on. The initial time period where these quantities evolve 
most rapidly (say t < 20000) corresponds to the first two stages of the dynamics described above, in which the 
modulational instability creates an array of soliton-like structures which then coalesce into a single coherent soliton. 
After the coalescence has ended, the kinetic (potential) energy increases (decreases) very slowly to its saturation value. 
In the process, fluctuations develop on finer and finer spatial scales, which accounts for the gradual increase of kinetic 
energy, while the surviving soliton slowly absorbs mass from the background fluctuations, thereby increasing the 
magnitude of the contribution to the potential energy from the coherent structure. In the long-time limit, therefore, 
the soliton accounts for the vast majority of the potential energy, while the fluctuations make a substantial contribution 
to the kinetic energy. 

The mean-filed statistical theory provides a prediction for the expected value of the kinetic energy Kn in statistical 
equilibrium for a given number of modes n. This is {Kn) = (('!/''"■*)) + — H*, which follows directly upon 
multiplying eqn. (23) by fc| and summing over j. The first term in this expression for (Kn) is the contribution to 
the mean kinetic energy from the coherent soliton structure which minimizes the Hamiltonian Hn subject to the 
particle number constraint iV„ ~ N'^. The second term in (A'„) is the contribution to the expectation of the kinetic 
energy from the fluctuations. H* is the minimum value of Hn given the particle number constraint. As n — > cx3, we 
see that (Kn) converges to K{ijj°°) + — H*, where is the minimizer of the Hamiltonian H given the particle 
number constraint N — for continuous NLS system on the interval [0, L], and H* = H{ip°°). Approximating 
K{ip°°) and H{ip°°) by K{(p) and H{ip), where (p is the solitary wave on the real line whose particle number is iV°, we 
obtain for the setting considered in Figure (^) the large n estimates Kn{{'4>^'^^)) ~ 9.2, — iJ* w 22.4, and therefore, 
{Kn) ~ 31.6. Also, according to the statistical theory, the expected value (0„) of the potential energy in statistical 
equilibrium should converge as n ^ oo to Approximating this by 6((^), with (f> as above, we have the estimate 

(0„) w —37.1, which we expect to be accurate for sufficiently large n. We see that the kinetic (potential) energy of 
the numerical solution is bounded above (below) by the estimate based on the statistical theory, but as expected, 
the solution does not attain the theoretically predicted value for a finite number of modes. This is because, for the 
spectrally truncated system, a finite amount of the particle number and the potential energy integrals are actually 
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contained in the small-scale fluctuations (according to the statistical theory, the contribution of the fluctuations to 
these quantities should be 0(l/n), where n is the number of spectral modes -this follows from ( p3|) IT^]). It may 
be checked that the spatial resolution is improved (i.e., when the number of modes n is increased, while the length 
L of the spatial interval, and the values and of the Hamiltonian and the particle number are held fixed), 
the contributions of the fluctuations to the particle number and the potential energy decrease, and the saturation 
values of the kinetic and potential energy attained in the numerical simulations come closer to the predicted statistical 
equilibrium averages of these quantities (see the inset in Fig. (||), which shows that the saturation value of kinetic 
energy increases towards the predicted statistical equilibrium value of modes n in the numerical simulation increases) . 
We expect that the contributions of the fluctuations to the particle number and the potential energy should vanish 
entirely as n ^ oo for fixed L, and A^°, and that the predicted statistical equilibrium values for the mean kinetic 
energy and potential energy should be approached very closely by the numerical solution in the long-time limit when 
the number of modes in the simulation is sufficiently large. 
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FIG. 4. Time evolution of the kinetic (upper curve) and the potential (middle curve) energies. The kinetic energy is 
increasing and consequently the potential energy is decreasing, in accord with the statistical theory developed above. The lower 
line indicates the potential energy of the solitary wave that contains all the particles of the system. The curves are obtained 
from an ensemble average over 16 initial conditions for n = 512. The length of the system is L = 128, and the (conserved) 
values of the particle number and the Hamiltonian are, respectively, A'^" — 20.48 and = —5.46. 

Figures (|l|) and (^) clearly illustrate that for a given (large) number of modes n, the dynamics converges when 
t ^ oo to a state consisting of a large-scale coherent soliton, which accounts for all but a small fraction of the particle 
number and the potential energy integrals, coupled with small-scale radiation, or fluctuations, which account for the 
kinetic energy that is not contained in the coherent structure. Formula (|2^) suggests, in fact, that in the long-time 
limit, the coherent structure and the background radiation exist in balance (or in statistical equilibrium) with each 
other, through the equipartition of kinetic energy of the fluctuations. In Figure |^), we display the particle number 
spectral density IV'fePi where i/'fe is the Fourier transform of the field ijj, as a. function of the wave number k for a 
long time run. To obtain this spectrum, we have performed both an ensemble average over 16 initial conditions, 
and a time average over the final 1000 time units for each run. For comparison, we have displayed in this figure the 
spectrum of the solitary wave ( p5| ) whose particle number is equal to conserved value of the particle number for the 
simulation. Observe that there is both a qualitative and quantitative agreement between the spectrum of this solitary 
wave solution and the small wavenumber portion of the spectrum arising from the numerical simulations. This is in 
accord with the statistical equilibrium theory, which predicts that the coherent structure should coincide with this 
solitary wave (in the limit n — )■ oo). For larger wavenumbers, the spectrum of the numerical solution is dominated 
by the small scale fluctuations. We have indicated on the graph the large wavenumber spectrum pre dicted by the 
statistical theory. This prediction comes from the second expression on the right hand side of eqn. (|2^), except that 
we have approximated the minimum value ff* of the Hamiltonian for the spectrally truncated system with n modes 
by the Hamiltonian H* of the above-mentioned solitary wave solution for the continuum system. Not only is there 
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a good qualitative agreement with the predicted equipartition of kinetic energy amongst the small-scale fluctuations 
(i.e., the fc^^ slope), but there is also an excellent quantitative agreement between the numerical results and the 
formula (|2^) for large k. Let us mention that the long-time spectrum obtained from a single simulation starting from 
a given initial condition, and without time averaging, though similar to the spectrum displayed in Figure (|^), is much 
noisier. 
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FIG. 5. Particle number spectral density [i/'fel^ as a function of k for t = 1.1 x lO'' unit time (upper curve). The lower curve 
(smooth one) is the particle number spectral density for the solitary wave that contains all the particles of the system. The 
straight line drawn for large k corresponds to the statistical prediction (|2^) for the spectral density for large wavenumbers. The 
numerical simulation has been performed with n = 512, dx = 0.25, Af" = 20.48 and H° = -5.46. 

As we have mentioned above, the numerical spectrum shown in Figure (H) arises from an ensemble average over 
long time and over different initial conditions (with the same values of the particle number and the Hamiltonian). 
Now, under the assumption that the dynamics is ergodic, such an average should coincide with an average with 
respect to the microcanonical ensemble for the spectrally truncated NLS system |l^. Since it can be shown that 
the the mean-field statistical ensembles p'-"-' constructed above concentrate on the microcanonical ensemble in the 
continuum limit n ^ oo (see Theorem 3 of reference [|l^), it should be that averages with respect to p*^"^ for large 
n agree with the ensemble average of the numerical simulations over initial conditions and time, assuming ergodicity 
of the dynamics. While we have not shown that the dynamics is ergodic, we have, in fact, demonstrated what we 
believe to be a convincing agreement between the predictions of the mean-field ensembles p*-"^ and the results of 
direct numerical simulations. 

We have also monitored the time evolution of the quantities 
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(26) 



for m a positive integer. For the periodic boundary conditions considered here, the index j ranges from —n/2 to 
n/2 and the wavenumber kj is given by kj = l-KjjL. Note that S\ is the kinetic energy. In general, Sm{'>P) is the 
squared L2 norm of the m-th derivative of the field "0. The growth of Sm in time is an indicator of the development of 
fluctuations of the fleld on fine spatial scales. In addition, we may consider that Sm gives an estimate of the evolution 
of the typical wave number K{t) of the fluctuations since, roughly speaking, we can estimate Sm K{ty^™~^\ 

The mean-fleld statistical theory provides the following prediction for the expectation of Sm in statistical equilibrium 
for a given number of modes n: 
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(27) 



where we have used eqn. (p3|). The first term is the contribution to {Sm) from the mean field (the coherent structure), 
and the second term is the contribution from the fluctuations. Note that for a finite number of modes n, {Sm) is finite 
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for each m, but only (Si), which is the mean of the kinetic energy, remains finite in the continuum hmit n — > oo. 
The divergence of (Sm) for m > 2 comes from the second expression on the right hand side of ( ^7| ) (i.e., from the 
fluctuations), which is of the order ^^(m-i) as n oo. For example, when m = 2 this expression is found to be 
7t^{H'^ — H*){n^ + 3n + 2)/(3i^), and we have the following formula for {S2) for a given number of modes n in the 
spectrally truncated NLS system: 



(82) 



n/2 



Tr\H°-H*)in^ + 3n + 2) 



3L2 



(28) 



Based on the considerations of the previous paragraph, we expect that the numerical simulations for given number 
of modes n will reveal that the quantities Sm are bounded and saturate when i — > 00, but that the larger the number 
of modes n, the larger the saturation value of Sm (at least for m > 2). Figure (^ shows the evolution in time of ^2 
for different values of n (with the same L, and iJ"). We observe that saturation does indeed occur for a finite 
number of modes. Also, as n increases, the saturation value increases, as does the time required to reach saturation. 
By approximating the sum in eqn. ( p8[ ) by \4'xx\'^ dx and approximating H* by H{(j)), where (f> is the solitary 
wave on the whole real line whose particle number is equal to the conserved particle number for the simulations 
treated in Figure (^), we obtain the following estimates: (52) « 27.1,45.5,97.6 and 170.3 for n = 48,64,96, and 
128, respectively. Note that these estimates for (^2) agree closely with the observed saturation values of 6*2 in the 
numerical simulations for n = 48 and n — 64. For n = 96, saturation has not quite yet been reached, but the value 
of 5*2 at the final time t = 3 x 10^ is still reasonably close to the theoretical estimate of 97.6. For n — 128, ^2 is still 
growing considerably at the final time of the simulation, and so we can not make comparisons with the statistical 
prediction for (^2) at this point. The inset in Figure (^) shows the evolution of the kinetic energy 6*1 as a function of 
time for n = 48, 64, 96, and 128. We see that the kinetic energy saturates nearly at the same rate for all of the values 
of n considered here. Clearly, S'l remains bounded as n increases. As discussed above, (^i), the statistical equilibrium 
value of of the mean kinetic energy, converges as n — > 00 to K{iIj°°) + iJ" — H{ip°°), where "0°° is the solitary wave 
that minimizes the Hamiltonian H for the given particle number N'^ for the NLS system on the interval [0, L]. Once 
again, approximating 'ip°° by the solitary wave (f) on the whole real line that minimizes H given the particle number 
constraint N = N^, we may estimate the limiting value (S'l) by K{(l)) + i?" — H{(f>), which is « 7.3 for the value 
N*^ — 9.6 considered in Figure (^. This estimate provides an upper bound on the saturation values of S'l observed in 
the simulations, and as the number of modes in the simulations increases, Si saturates closer to this approximation 
of the statistical equilibrium value. 
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Time xio' 
FIG. 6. S2 as a function of time for n = 48, 64, 96 and 128 (lower to upper curve). The simulations are all performed for a 

box of size L = 38.4, for iV" = 9.6 and H° = -3.2. The curves are obtained from an ensemble average over 16 runs for each 

n. Saturation is reached for n = 48 and n = 64, while it is almost obtained for n = 96. As n is increased, the time required to 

reach saturation increases, and the saturation value also increases. The inset shows the kinetic energy Si as a function of time 

for the same values of n. As opposed to ^2, the saturation of Si seems to be occurring at about the same rate for each n. The 

saturation value of Si increases slightly as n increases, but it remains bounded above by the statistical equilibrium value (Si). 
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When the spatial resolution of the numerical simulations is improved (i.e, when n is increased with L fixed), the 
functions Sm are typically seen to exhibit power law growth in time before reaching saturation (see Figure |^). 
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FIG. 7. Sm. as a function of time for m = 2,3 and 4 (lower to upper curve). The growth of these quantities is indicative 
of the development of fluctuations on fine spatial scales as time increases. The Sm have been calculated for 512 modes with 
N° = 12.8, H° = -4.26 and dx = 0.1, with an ensemble average over 16 runs. They have been rescaled in order to display 
them on the same graph. 



Indeed, we observe for m — 2,3 and 4 that Sm obeys the following power law dynamics: 



(29) 



with v = 0.25 ± 0.01. This behavior is observed for t large enough that the coalescence process has ended. It 
corresponds, therefore, to the regime where the kinetic energy has essentially reached saturation. Remarkably, the 
observed dynamical exponent ly is in good agreement with the prediction of Pomeau fl^ , which estimates the evolution 
of the typical wave number K{t) of the fluctuations as time increases. The estimate comes from a dimensional analysis 
of the weak turbulence equation deduced from equation (|23). Describing the fluctuation field as: 



the relation between the energy w and the wave number k is called the spectrum of excitations (we refer the reader 
to |l^ for details). In jl^, it has been shown that if this wave number K{t) is in the range where the spectrum of 
excitations obeys u){k) ~ fc^, which means that the fluctuations behave essentially like free particles, then, assuming 
that there is two-wave resonance in the weak turbulence approximation, it follows that K(t) ^ (ef)^^^. Here, e is 
the spatial energy density of the fluctuations (so e ^ (i?" — H*)/L). The analysis in [Q was carried out for cubic 
defocusing NLS in two spatial dimensions, and does not immediately go over to the case under consideration here. 
In fact, strictly speaking, in the regime to [k) = k'^, the resonance of two waves cannot hold in one spatial dimension. 
However, we conjecture that, due to the interactions with the large-scale coherent structure, the resonance may in 
fact be meaningful in the present setting, and therefore, we believe that a dimensional weak turbulence analysis along 
the lines of that developed in may be relevant. We hope to explore this possibility in the future. Interestingly, 
for the NLS system, the approximation oj{k) = k'^ is usually valid in the limit fc ^ 1. In the numerical simulations, 
the finite number of modes provides an ultraviolet cutoff since the largest wave number of the system is 

TT TT 

dx L 

We remark that we have been able to recognize the power-law growth in time of the quantities Sm only for the smallest 
dx we have considered in our simulations. For larger dx the free particle regime might not be realized, and it is not 
surprising in this case that the power law behavior is not observed. 
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The previous considerations allow us to attach a more precise meaning to what we have been referring to as the 
transitory dynamics and the equilibrium state for the spectrally truncated NLS system. For sufficiently small grid 
sizes dx, one may consider that Sm grows according to the power law (|2^) until the time tn at which the typical wave 
number K{t) reaches the largest available wave number kmax- This time i„ (which appears to scale as n'* for fixed 
L, TV" and H^) defines a crossover between the transitory regime in which the solution evolves towards the spectrum 
( P3| ) , and the statistical equilibrium regime where the system investigates its phase-space according to the probability 
density p^"\ Notice that in the contimmm limit n ^ oo, tn diverges to infinity, so that continuous NLS system can 
not reach statistical equilibrium in finite time. Such conjectures are supported by the investigation of the dynamics 
of the particle number spectrum during the intermediate time regime after which the coarsening process has ended, 
but before the final statistical equilibrium state has been reached. Note that the statistical equilibrium model does 
not provide predictions about the time evolution of the spectrum, because it is strictly an equilibrium theory. In fact, 
based on the statistical theory alone, nothing can be said about path which the system follows from the statistically 
unlikely initial condition to the final statistical equilibrium state. Figure (^) displays the particle number spectrum 
at t = 5 • 10^ unit time for a spatial resolution of da; = 0.1. 




wave number k 



FIG. 8. The particle number spectral density for n — 512 and dx = 0.1 (thus L — 51.2) at unit time t = 5 - 10^. The coherent 
soliton structure already accounts for almost the entire number of particles of the system, but the system has not yet reached 
statistical equilibrium. The initial noise level is still present for large wavenumbers {k > 20), while at smaller wavenumbers, 
one can recognize both the soliton-like structure and a fluctuation spectrum following approximately a law. The spectrum 
has been obtained by an ensemble average over 16 initial conditions and a time average over the final 10 unit times. 



This figure illustrates that the system investigates smaller and smaller scales as time increases. Indeed, the particle 
number spectrum for the modes A: > 20 is still at the level of the initial noise. Thus, the smallest scales available 
to the system have yet to be excited at the time t = 5 ■ 10^. For larger scales, however, one can recognize both the 
spectrum corresponding to the coherent soliton structure, and the fluctuation spectrum which appears to follow, at 
least approximately, the equilibrium law IV'fcP k~^. This suggests the following scenario for the spectrally truncated 
NLS dynamics: as time increases, smaller and smaller scales are explored, until eventually all available modes are 
excited. However, at any given time after the coalescence process has ended, the system may be considered as being 
in statitistical equilibrium over all the modes that have been excited up to that point. Defining kmaxit) to be the 
largest wave number that the system has reached up to time t, then we obtain from (|2^) that kmax{i) ^ (^O^^^i if 
kmax is large enough. Now if we denote by n{t) the number of modes that have been excited up to time then based 
on our previous arguments, we have that kmax{t) — TTn{t)/L. But, using kmax{t) ~ (e^^^^'i obtain the following 
estimate for the spectrum of the fluctuations at time t: 

ffO _ TT* 3/4 

for < kyna.x{t). The total particle number spectrum at time t, of course, has to be taken as the sum of the spectrum 
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corresponding to the large-scale coherent structure (which decreases exponentially for large k) and the spectrum of 
the fluctuations. 

We emphasize that the derivation of the equation ( ^ ) describing the particle number spectrum at an intermediate 
time t crucially depends on the assumption that the system evolves in such a way that it is nearly in statistical 
equilibrium over all the modes that have been investigated up to that time. The Figure has motivated us to 
make this assumption, but clearly further numerical investigations should be carried out in order to test the validity 
of this hypothesis, as well as the accuracy of the formula (|30|). Nevertheless, we find it quite interesting that the 
fluctuation spectrum (^ ) agrees with the prediction in which was derived from a dimensional analysis of the 
weak turbulence equations for the NLS system. 

V. DISCUSSION-CONCLUSION- ACKNOWLEDGEMENTS. 

The primary purpose of the present work has been to test the predictions of a mean-field statistical model of self- 
organization in a generic class of nonintegrable focusing NLS equations defined by eqn. (|l|). This statistical theory, 
which has been summarized above, was originally developed and analyzed in jist . In fact, we have demonstrated a 
remarkable agreement between the predictions of the statistical theory and the results of direct numerical simulations 
of the NLS system. There is a strong qualitative and quantitative agreement between the mean field predicted by the 
statistical theory and the large-scale coherent structure observed in the long-time numerical simulations. In addition, 
the statistical model accurately predicts the the long-time spectrum of the numerical solution of the NLS system. 
The main conclusions we have reached are 1) The coherent structure that emerges in the asymptotic time limit is the 
solitary wave that minimizes the system Hamiltonian subject to the particle number constraint N = N'^, where 
is the given (conserved) value of N, and 2) The difference between the conserved Hamiltonian and the Hamiltonian 
of the coherent state resides in Gaussian fluctuations equipartitioned over wavenumbers. 

While the statistical model we have developed is an equilibrium theory, and, strictly speaking, only provides 
predictions concerning the long-time statistical properties of the NLS system, we have combined this theory with 
insight gained from numerical simulations to paint a picture of the nature of the dynamics leading up to the statistical 
equilibrium state. Specifically, the simulations (and, in particular, the results shown in Fig. ^) indicate that evolution 
is such that, at a given time after the coarsening process has ended, the system is nearly in statistical equilibrium 
over all the modes that have been excited by that time. From this observation the fact that the quantities S',„ defined 
in (p6|) are seen to exhibit the power law growth in time according to eqn. (|2^), we have arrived at the prediction 
( ^ ) for the time-dependence of the spectrum of the fluctuations. As we have mentioned above, results such as ( |29| ) 
and ( |30| ) have previously been derived by Pomeau [ p^ from weak turbulence arguments, but for the defocusing cubic 
NLS equation in a bounded two-dimensional spatial domain. We believe that it would be an interesting exercise to 
check whether these formulas can be derived directly from a weak turbulence analysis in the present context -that is, 
for ID nonintegrable focusing NLS equations in the absence of collapse. 

We would like to point out certain analogies between the dynamics of the NLS systems we have considered here, and 
the dynamics of a turbulent 2D Navier-Stokes fluid. A prominent feature of large Reynolds number 2D Navier-Stokes 
turbulence is the formation of quasisteady coherent vortex structures [0-1). Starting from generic initial conditions, 
the evolution of the fluid is characterized by the formation of a collection of large-scale vortices, and the subsequent 
merger or coalescence of like-signed vortices ^J. The large-scale soliton structures in (focusing, nonintegrable) NLS 
play a role similar to that of the vortices in 2D Navier-Stokes turbulence. Indeed, we have observed in our numerical 
simulations of NLS the formation of an array of soliton-like structures which eventually coalesce into a single persistent 
soliton of large amplitude. Another characteristic feature of turbulence in two-dimensions is the presence of a dual 
cascade Q. There is a direct cascade of enstrophy to small scales and an inverse cascade of energy to large scales. As 
pointed out long ago by Kraichnan [Q, the existence of the inverse cascade of energy is indicative of the formation 
of a large-scale structure in the system. In NLS, there is also a dual cascade. Indeed, our numerical simulations, 
which correspond to injecting as initial conditions particle number and energy at a given scale h associated with the 
modulational instability, have revealed that their is a direct transfer of kinetic energy to spatial scales smaller than 
li, while the particle number is transferred to large scales. While the ID NLS equation is much simpler system to 
investigate, both analytically and numerically, than a turbulent 2D fluid system, we believe that the understanding 
of the coalescence and transfer processes in this generic model of nonlinear wave turbulence might provide important 
insight into the nature of turbulent systems in general. 

It is a pleasure to thank Robert Almgren, Shiyi Chen, Leo Kadanoff, Yves Pomeau, Bruce Turkington, and Scott 
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